Effect of thermal fluctuations in topological p-wave superconductors 
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We study the effect of thermal fluctuations on the topological stability of chiral p-wave supercon- 
ductors. We consider two models of superconductors: spinless and spinful with a focus on topolog- 
ical properties and Majorana zero-energy modes. We show that proliferation of vortex-antivortex 
pairs above the Kosterlitz-Thouless temperature Tkt drives the transition from a thermal Quantum 
Hall insulator to a thermal metal/insulator, and dramatically modifies the ground-state degeneracy 
splitting. Therefore, in order to utilize 2D chiral p-wave superconductors for topological quantum 
computing, the temperature should be much smaller than Tkt- Within the spinful chiral p-wave 
model, we also investigate the interplay between half-quantum vortices carrying Majorana zero- 
energy modes and full-quantum vortices having trivial topological charge, and discuss topological 
properties of half-quantum vortices in the background of proliferating full-quantum vortices. 
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I. INTRODUCTION 

Topological phases of matter have been subject of in- 
tense physics research in the last decade. 1 In addition to 
interest from the fundamental physics perspective, these 
states of matter can also be used for topological quantum 
computation, which is predicted to have an exceptional 
fault-tolerance by virtue of encoding and manipulating 
information in non-local degrees of freedom of topolog- 
ically ordered systems. 3 Candidate physical systems in- 
clude Fractional Quantum Hall states 1 and topological 
superconductors. • In all these systems the topological 
degrees of freedom coexist with non-topological ones. It 
is important to understand their interplay because it of- 
ten determines the stability of the topological phase. 

In this Letter we focus on 2D topological p-wave super- 
conductors and study their robustness against thermal 
fluctuations. Specifically, we investigate the topological 
degeneracy in these systems in the presence of thermally- 
generated topological defects (vortices). Without vor- 
tices, the stability condition for the topological supercon- 
ducting phase is set by the quasiparticle energy gap A, 
i.e., T <C T c ~ A. We will show that vortex-antivortex 
proliferation provides a more stringent temperature re- 
quirement. 

We first focus on a spinless model, where we model 
thermal fluctuations with a classical XY model. Increas- 
ing the temperature above the Kosterlitz-Thouless (KT) 
transition point but still well below the local quasipar- 
ticle gap (Tkt < T -C T c ), vortices start to prolifer- 
ate and eventually destroy the topological spinless super- 
conducting phase by driving the system into a thermal 
metal or non-topological insulator phase. The degen- 
eracy splitting in the low-temperature (T < Tkt) and 
high-temperature (T > Tkt) phases changes from an ex- 
ponential to a power-law scaling in the system size. 

A spinful model allows for both half-quantum (HQV) 
and full-quantum (FQV) vortices. Only the former carry 
robust Majorana zero-energy modes. Thus, one can con- 



sider the interesting situation where the superconduct- 
ing phase is disordered due to the presence of FQVs and 
study the splitting of a degeneracy due to HQVs embed- 
ded in the system at large enough separation R. Naively, 
one might expect that the splitting would not be affected 
by the proliferation of FQVs since the splitting energy 
is governed by the local quasiparticle gap which is only 
weakly affected by thermal fluctuations. However, we 
show that the situation is much more intricate and re- 
quires a deeper understanding of the interplay between 
topological and non-topological degrees of freedom. 

Disorder in superconductors in symmetry class can 
drive a transition from the thermal Quantum Hall (TQH) 
phase to either a thermal metal (TM) phase or a topolog- 
ically trivial thermal insulator (TTI). 5,9-13 These phases 
have been primarily studied within network models 1 1-1 3 
and very recently in certain microscopic models. 1 ,15 We 
find that all these phases appear in our microscopic model 
with thermally generated disorder. The TTI is an Ander- 
son insulator with a non-zero density of states at E = 0. 
We demonstrate that it does not realize the aforemen- 
tioned scenario where the splitting is governed by the lo- 
cal quasiparticle gap if the phase is disordered. Instead, 
the presence of additional emergent zero-energy modes 
changes the splitting of the topological degeneracy from 
exponential to power-law. 



II. MODEL 

We set up our numerical problem in two steps. First, 
we model thermal fluctuations of the superconducting 
phase using a classical XY Hamiltonian. We then make 
an adiabatic approximation assuming that the vortex 
dynamics are slow compared to the quasiparticle one 
which appears to be quite reasonable since Abrikosov 
vortices are macroscopic objects and have large effec- 
tive mass. Under these conditions, quasiparticles are 
moving in a static background of different vortex con- 
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FIG. 1. Schematic phase diagram of a spinful p-wave super- 
conductor as a function of temperature and in-plane magnetic 
field. The phase diagram of the spinless case is recovered for 
By — 0. Above (By) c , the superconducting order parameter 
vanishes in a self-consistent calculation. 

figurations. Secondly, we diagonalize the Bogoliubov- 
de-Gennes Hamiltonian for each disorder realization and 
compute the quasiparticle energy spectrum, the density 
of states (DoS) and inverse participation ratio (IPR). 

We consider a model for a p + ip superconductor on a 
torus of L x L sites defined as 

H = E ^ a '4a^+^f4J^+h.c) (1) 

(i,j),a,cr' 

We first study a spinless model by choosing t aa = 
—tS^a-i, Af? — AoXijdijSa-^i , where \x is the chemical 

potential. This corresponds to a d-vector characterizing 
spin-triplet pairing to be aligned along ai-axis. With- 
out spin-mixing perturbations we can equivalently study 
a spin-polarized system. Xij is a chirality factor that im- 
plements p + ip pairing and is ±1 for j = i ± x, and ±i 
for j = i ± y. 9ij is a phase variable to be discussed be- 
low, and Ao is chosen to be a constant. We will discuss 
the self-consistency condition for Ao below. We solve 
the corresponding BdG equation numerically to obtain 
eigenvalues E n and eigenstates (u n ,v n ) T . 

In the homogeneous case 9ij = 1, the dispersion of the 
spinless Hamiltonian is 

e(k) = — fi— 2t(cos(k x ) + cos(k y )) (2) 

E{k) =v /e(fc) 2 + Aotsin 2 ^) + sin 2 (*;„)). (3) 

The dispersion has gapless points for fi = — At, 0, At. The 
homogeneous system is in a topological phase if Ao > 
and -At < fi < or < fi < At. 

Vortices in a spinless p-wave superconductor bind Ma- 
jorana zero-energy modes. These localized quasiparti- 
cles are described by a self-conjugate operator 7 = 7^. 
The ground state degeneracy and Majorana quasiparti- 
cles lead to non-Abelian braiding statistics in these many- 
particle systems. 18-21 Depending on parameters, vortices 



may carry a large number of localized states below the 
bulk gap. To reduce the required computational effort, 
we typically use A ~ fi. For such a choice of parameters, 
the coherence length \ becomes comparable to the lat- 
tice spacing and no midgap states except the zero-energy 
state are present. The effect of such midgap states has 
been treated in Ref. 22 and will be discussed at the end 
of our paper. 

The situation becomes more subtle when many vortices 
are present and localized zero-energy modes hybridize 
leading to a ground-state degeneracy splitting. 23-26 At 
large vortex separation R 3> £, Majorana modes acquire 
an exponentially small energy splitting and the ground- 
state degeneracy at small vortex density is preserved up 
to exponentially small corrections SE, which for p-wave 
superconductors read ' 

.„ [SAT? ( A 2 \ 1/4 Y(kR) ( R\ 

2 

Y(kR) =cos(fci? + a) - -sin(kR + a) (5) 
A 

2(1 + A 2 ) 1 / 4 

H 

+ A 

with A = fc£, 2a — arctan A, k = \j2rapL — Aq/u|,, the 
Fermi velocity vf and superconducting coherence length 
£. Thus, the effective low-energy model for a multi- vortex 
configuration reads H = ^^ijlilji where jj is a 

self-conjugate (Majorana) operator representing a zero- 
energy state in j-th vortex. Given that in realistic sys- 
tems kRij 3> 1, SEij is a rapidly oscillating function. 

We now study how this ground-state degeneracy is 
modified by vortex-antivortex proliferation above the KT 
transition. We consider a situation where thermal fluctu- 
ations affect only the phase of the order parameter while 
the magnitude remains approximately constant. Below, 
we will give evidence from a self-consistent calculation 
that such a regime can be obtained. The fluctuations of 
the phase can be modeled by a classical XY Hamiltonian 
for the phases 9ij, 

H = -J^2 cos(arg% - arg6> ifc ), (6) 

where J is related to the superfluid stiffness. A key prop- 
erty of this model is that below the KT temperature 
(Tkt = 0.89 J in the infinite system 28 ' 2!) )vortices and an- 
tivortices are bound in pairs by a logarithmic attraction. 
Above the transition, they unbind and proliferate. The 
Monte Carlo sampling is performed using a standard clus- 
ter update method. ' 

To study the effect of thermal fluctuations on 
the topological degeneracy, we introduce a fixed vor- 
tex/antivortex pair in the system by adding a non- 
fluctuation phase factor to A,^ and study the energy 
splitting in the presence of the background defects. The 
fixed vortex-antivortex pair is implemented by introduc- 
ing an additional phase factor, 

% = exp (ifcj) = exp (ififj) exp (-i<f>ij) ■ (7) 
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FIG. 2. (color online) Self-consistent calculation for A for 
lattice of size L = 16 and [7 = 5. The Kosterlitz-Thouless 
transition takes place at T/ J — 0.89 and leads to the first drop 
in the gap A, which renormalizes it quantitatively while the 
system remains in a superconducting phase. For even higher 
temperatures, the gap vanishes and the superconductivity is 
destroyed. 



FIG. 3. (color online) Density of states D(E) in the low- 
temperature phase (T = 0.1,0.3,0.5,0.7,0.9, from right to 
left) for L — 64. The DoS is strongly suppressed for suffi- 
ciently low temperatures and only shows a peak close to zero 
energy corresponding to a slightly renormalized value of the 
energy of the T = case, and an increase as the transition 
is approached. Simulations were performed for Ao = t/2 and 
/i = —t. 



III. RESULTS 



Here {4>fj) are the polar angles that the bond ij 
has with the vortex (antivortex) located at position A 
(£>). When applying this to a torus mapped to a lattice 
with periodic boundary conditions, special care has to 
be taken that the order parameter is smooth around the 
boundary. 



To obtain a simpler description on the torus, we per- 
form a gauge transformation after which the vortices 
are implemented only by a tt phase shift in both the 
hopping and the pairing terms across a particular line 
(branch cut) connecting the two vortices. To this end, 
we introduce a gauge field <pi on the sites such that 
4>ij = {4>i + 0j)/2 and perform the gauge transforma- 
tion c, — > Ci exp (— i(j>i). For a single vortex, the field 4>i 
would wind from to n around the vortex, whereas 4>ij 
would wind from to 2tt; therefore exp{i<j)ij) is smooth, 
while exp(— ifa) has a jump from —1 to 1 at a branch cut. 
It follows that everywhere except along this branch cut, 
exp(— exp(—i(j)j) exp(—i(j)ij) = 1 and the phase is re- 
moved from the anomalous hopping; along this branch 
cut, a phase 7T remains and hopping terms with i and 
j on different sides of the branch cut pick up a minus 
sign. For the situation of two vortices, this can be gener- 
alized and one finds that the branch cut turns into a line 
connecting the two vortices. Applying the gauge trans- 
formation in the normal hopping, one finds that the same 
phase factor across the branch cut is introduced. 



A. Self-consistent calculation 

Owing to the mean-field approximation, the BdG 
equations have to be augmented with self-consistency 
conditions, which read 

Aij = U(c t c 3 ) (8) 
= UJ2 <(*HWtanh(g^ (9) 
A = (|Ay|) (spatial and MC average). (10) 

We have performed self-consistent calculations for a 
spinless p-wave superconductor without fixed vortices 
and obtained the order parameter A = \J (ciCj) 2 as a 
function of temperature. Our results are shown in Fig. 2 
for different values of J, which controls the relative tem- 
perature scales of the Kosterlitz-Thouless transition and 
the mean-field transition where superconductivity is de- 
stroyed. Our data show that for sufficiently small J, these 
two transitions are well-separated and there is an inter- 
mediate regime where the phase of the order parameter 
is disordered, but its magnitude remains finite at a value 
20-30 % below the zero-temperature result. Since we are 
interested in qualitative results only, such a small quali- 
tative change is irrelevant and we do not perform a self- 
consistent calculation but fix a value of A independent 
of temperature. 
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B. Low-T phase 



Diagonalizing the BdG equation for each configuration 



of dij, we compute the DoS 



D(E) 




(ii) 



where (-)t indicates the Monte Carlo average at tem- 
perature T, and N denotes the number of states. We 
generally average over at least 10,000 configurations and 
obtain error bars with a standard Jackknife analysis. As 
shown in Fig. 3, the DoS at zero temperature shows a 
sharp peak at the energy splitting set by the system size 
for the fixed vortices, and a continuum of states above 
the bulk gap A . At low temperatures T <C 7kt> both 
features are broadened but the energy splitting of the 
Majorana modes remains exponential and the DoS is sup- 
pressed between this scale and the bulk gap. 

To further elucidate the fate of the ground-state degen- 
eracy, we study the energy splitting between fixed vor- 
tices by fitting it to (cf. Eqn. (4)) 



oh — — = exp 



(1 + c 2 cos(c 3 a; + c 4 )) , (12) 



where £, Ci, C2, C3 and C4 are fit parameters. Good fits 
are obtained for low temperatures, as shown in the in- 
set of Fig. 4. Our results for the correlation length are 
shown in the main panel of that figure. The correlation 
length depends only weakly on temperature as long as 
the system is well below the KT transition. At the tran- 
sition, a sharp jump in the correlation length indicates a 
fundamental change in the scaling behavior. 



C. High-T phase 

The sharp change is related to the delocalization tran- 
sition (i.e. appearance of a disorder-driven thermal metal 
phase) characteristic to class D superconductors. 9 ' 12 ' 13 ' 15 
The TM is characterized by delocalized states at E = 
and a logarithmic divergence of the DoS for low 
energies. 15,31 Furthermore, the oscillatory behavior of the 
DoS in the zero-dimensional limit is consistent with the 
random matrix theory predictions for class D. In Fig. 5, 
the DoS for a spinless p-wave superconductor well above 
the KT transition is shown along with a fit to the random 
matrix theory result s for symmetry class D 



D(E)~ 1 - 



sm(2Tr~fEL 2 
2ttEL 2 



(13) 



Using a single-parameter fit, we obtain excellent agree- 
ment with our theoretical expectations for the TM phase: 
i) At the lowest energies, the DoS follows random matrix 
theory predictions, ii) For higher energies (but still well 
below the bulk gap A), a logarithmic divergence is ob- 
served. This clearly establishes that there is a TM phase 
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FIG. 4. (color online) Scaling of the energy splitting with 
vortex distance in a system of size L = 64. The main panel 
shows the dependence of the coherence length £ on the tem- 
perature. A KT transition takes place at T ~ 0.9. The inset 
shows the dependence of the splitting e on the distance with 
a fit to Eqn. (12). Simulations were performed for Ao = t/2 
and /J, = —t. 




FIG. 5. Density of states D(E) in the high-temperature phase 
(T = 1.5, A = t/2, n = -t). Fits are to Eqn. (13). The data 
has been rescaled such that D(0.03) = 1. 



above the KT transition, as shown along the B y 
in Fig. 1. 



line 



The coefficient 7 above is related to the effective band- 
width in the Majorana fermion hopping problem defined 
above and should therefore be related to 6E (4). In- 
deed, we numerically confirm that 7 ~ Aoexp(— Aq/C) 
for some constant C. The energy scale is largely inde- 
pendent of temperature as long as the temperature is 
sufficiently far away from the KT transition. 
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FIG. 6. (color online) Top panel: Density of states for spinful 
fermions. With magnetic field, the DoS has a constant value 
at E — 0, whereas without magnetic field it diverges due to 
the thermal metal. In the case with magnetic field, a peak 
appears at E = B v . Bottom panel: Constant term obtained 
from a fit to I ~ a + L b . Values of a — > and a/0 indicate 
extended and localized states, respectively. 
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D. Spinful model 

When adding perturbations that break SU(2) symme- 
try, HQVs still carry zero-energy modes, but FQVs do 
not. For example, a Zeeman splitting generated by an in- 
plane magnetic field B y ^ (ic^cn — icj^c^ will move 

the the lowest excitation energy supported by an FQV to 
non-zero Eq = B y . Indeed for each pair of wavefunctions 
of the spinless model at energies ±2?, there are four wave- 
functions at energies ±E ± B y for this special choice of 
field direction. Thus, if the system has a band of delocal- 
ized states near E = 0, the system will remain in a TM 
phase for B y smaller than the width of this band and 
will transition to an insulating phase once B y is larger 
than the width of this band. To determine this width, in 
Fig. 6 we show the DoS and the localization properties 
of the states, which we characterize by the inverse par- 
ticipation ratio (IPR) 1(E) defined to be the average of 
the fourth moment of a wavefunction of energy E. The 
inverse participation ratio can be calculated in our setup 
using 



((u n ) 2 + ( V „) 2 ) 2 



S(E - E n 



(14) 



For finite systems, the IPR must be calculated by aver- 
aging over states in a finite range of energies centered 
around E by broadening the 5 function. 

For extended states, the dominant scaling of the IPR is 
expected to be 1(E) ~ L 2_!/ , where v is a non- universal 



FIG. 7. (color online) Fit of the IPR to I ~ a + L b for several 
parameters sets. The top panel shows B y = 0, T = 1.5; the 
center panel shows B y = 0, T = 3, and the bottom panel 
shows By = 0.25, T = 1.5. At T = 1.5, there is clearly only a 
very narrow band of delocalized states around E ~ B y . For 
T = 3, however, the band is broadened such that states at 
almost all energies appear delocalized for the system sizes we 
can access. 



correction to the exponent !J while the IPR is expected to 
approach a constant for localized states, limL-^-x, 1(E) > 
0. Therefore, in the thermal metal phase of the spinless 
model (or equivalently the spinful model at B y = 0), we 
expect 1(E) to scale with a power law at zero energy due 
to the presence of extended states. At non-zero energies, 
it is expected to approach a constant value for L — >• oo. 
In the spinful case at finite magnetic field, we expect 
delocalized states, and hence a power-law scaling of the 
IPR, at E ~ B y , while we expect localized states and 
a saturation of the IPR for other energies E < A. In 
all cases, we expect it to behave with a power law for 
energies higher than the local gap A. 

In the bottom panel of Fig. 6, our results obtained from 
the extrapolation of the IPR are summarized. These re- 
sults indicate the presence of delocalized states at energy 
E ~ By and E > A, as expected. A more detailed per- 
spective is provided in Figure 7, which shows several fits 
for the IPR for the spinful model with two different val- 
ues of the magnetic field, B y — and B y = 0.25, and two 
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FIG. 8. (color online) Main panel: correlation length £ for 
HQVs extracted from an analogous fit to Eqn. (12). The 
correlation length displays a clear jump at the KT transition. 
Inset: DoS for a system in the thermal insulator phase (B y = 
0.25, A = t, n = -2t, L = 40,48,56,64, T = 1.5) with 
(dashed lines) and without (solid lines) fixed HQVs. 
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FIG. 9. Contribution of two fixed half-quantum vortices to 
the density of states in the high-temperature phase of a spin- 
ful superconductor with finite magnetic field. The main panel 
shows the density of states (solid lines: with HQVs, dashed 
lines: without HQVs). The scaling of the characteristic en- 
ergy scale Ehqv, defined below Eqn. (15), is shown in the inset 
on a log-log scale. 



different temperatures T = 1.5 and T = 3. The constant 
terms extract from such fits for T = 1.5 are shown in 
Fig. 6. The top panel of Fig. 7 (B y = 0, T = 1.5) clearly 
shows the saturation of the IPR for energies < E < A, 
whereas for very small energies, such as E = 0.03, no 
clear sign of saturation is observed for the accessible sys- 
tem sizes. The middle panel shows the same situation 
for a higher temperature, T = 3. In this case, the IPR 
appears to follow a power law also for intermediate ener- 
gies such as E — 0.09. This is indicative of a delocalized 
band of finite width centered around E — 0, with the 
bandwidth growing as T is increased. Finally, the bot- 
tom panel shows the situation with finite magnetic field 
and temperature close to the KT transition, where no lo- 
calization is observed around E = B y , i.e. the IPR scales 
with a power law. 

We have checked that our result is robust against other 
perturbations, such as an additional magnetic field B x . 
Since our results are restricted to finite-size systems, it 
remains open whether these states are truly extended in 
the thermodynamic limit and how this can be connected 
to theoretical work. 



E. Half-quantum vortices 

We now consider the Majorana zero-energy modes car- 
ried by half-quantum vortices (HQVs) in the trivial ther- 
mal insulator (TTI) phase by studying the DoS with two 
fixed half-quantum vortices in the background of ther- 
mally fluctuating full quantum vortices. As shown in the 
inset of Fig. 8, the HQVs give an additional contribu- 
tion to the DoS at low energies. Studying the energy 
splitting as a function of the temperature, we find that 



the correlation length £ changes qualitatively at the KT 
transition, see Fig. 8. While our small system sizes do 
not let us determine whether the splitting is still expo- 
nential or becomes power law in this regime, the splitting 
energy for HQVs clearly changes dramatically above the 
KT transition. 

As an alternative approach to quantify the energy scale 
below which half-quantum vortices contribute, we study 
the integrated difference between the density of states 
with, D h (E), and without, D(E), half-quantum vortices: 

p(E)= [ E ' dE'(D h {E')-D{E')) (15) 
Jo 

and define Eh qv as the lowest energy such that 
Np(E hqv ) = 1, where N — 2L 2 . 

Fig. 9 shows the density of states with and without 
half-quantum vortices (cf. inset of Fig. 8). In the inset, 
the scaling of this quantity with system size is shown. 
A power-law scaling is clearly observed. This is consis- 
tent with the observation that the expectation value of 
the lowest energy (E )t behaves like a power-law both 
with and without HQVs. This shows that the topologi- 
cal degeneracy is destroyed by the presence of zero-energy 
states due to disorder, even though these states are lo- 
calized. 



IV. CONCLUSIONS 

We have studied the effect of thermal fluctuations on 
two-dimensional chiral p-wave superconductors in sym- 
metry class D. We have shown that thermally disorder- 
ing the superconducting phase drastically changes the 
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topological properties of these systems. We can ex- 
plain the underlying mechanism as the proliferation of 
vortex-antivortex-pairs, which carry low-energy excita- 
tions in their cores. Hybridization of these low-energy 
states gives rise to a thermal metal phase. Exploring the 
full phase diagram as a function of temperature and in- 
plane magnetic field, we find that in addition to the low- 
temperature topological superconductor and the thermal 
metal phase, a trivial thermal insulator (Anderson insula- 
tor) phase appears. We study the fate of the topological 
degeneracy in all these phases and find that it is well- 
defined only in the topological superconductor phase at 
T < Tkt- I n the thermal insulator phase, the splitting 
of the degeneracy due to half-quantum vortices changes 
dramatically due to vortex disorder. This result can 
be anticipated considering previous analytical evidence 
for a disorder-driven quantum phase transition in one- 
dimensional analogues of our system. 33 ' 3 

Throughout this paper, we have assumed A ~ ji in 
order to keep the coherence length on the order of a few 
lattice sites. In experimental systems, however, A is of- 
ten much smaller than ijl. In this case, one expects to find 
many subgap states localized on each vortex with energy 
splitting e s ~ A 2 //x. 22 While this regime is difficult to 
treat numerically, we argue based on our results for the 



high-T phase in magnetic field that there are two sce- 
narios for the high-T phase in the presence of many sub- 
gap states: if the hybridization scale for Majorana modes 
obeys SE < e s , the subgap states will contribute to the 
delocalized states in the thermal metal band which is 
therefore enhanced. If, on the other hand, e s SE, one 
would expect several peaks in the DoS centered around 
the subgap state energies £ s • n similar to the peak around 
the magnetic field. Here n is an integer corresponding to 
n-th energy level in the vortex core. In either case, ther- 
mal fluctuations destroy the topological phase and our 
qualitative conclusions regarding the disordered high-T 
phases remains valid. 
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